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Introduction 

Computational fluid dynamics (CFD) methods 
and advanced turbulence models are needed to pre- 
dict propulsion aerodynamic effects in transonic and 
supersonic free-stream conditions. Analytical results 
are frequently used to supplement the experimental 
data in critical design decisions. Accurate predic- 
tion of the pressure distribution and the skin friction 
coefficient is paramount to the design of propulsion 
systems. In the area of propulsion integration, ac- 
curate predictions of boundary layer structure, skin 
friction, and flow separation by CFD methods are 
critical. 

Two-equation turbulence models (ref. 1) offer sev- 
eral advantages over other approaches that compute 
practical flow problems. For example, algebraic mod- 
els lack turbulence history-dependent nonlocal effects 
(through the convection and viscous diffusion of the 
Reynolds stress models), effects which are known to 
be important in determining the turbulence struc- 
ture in complex flows. The numerical calculations 
that use the more advanced Reynolds stress mod- 
els (refs. 2 and 3) require the solution of trans- 
port equations for each component of the Reynolds 
stress tensor in addition to solution of the Navier- 
Stokes equations; this approach requires tremen- 
dous computational time for three-dimensional flow 
problems. The transport equations for second-order 
models require closure approximations for higher or- 
der turbulence correlations with uncertain physical 
foundations. Large-eddy simulations (ref. 4) con- 
stitute three-dimensional time-dependent computa- 
tions that require enormous computational time com- 
pared with traditional transport equation turbulence 
models. Moreover, the application of large-eddy sim- 
ulations to practical flows is often hindered by the dif- 
ficulties in modeling the turbulence near solid bound- 
aries and the problem of defiltering the results in 
complex geometries. 

Transport equations have been included in the 
standard, two-equation turbulence model of energy 
and dissipation rate (k-e). (See ref. 1.) The k-e 
equations can be applied to the near-wall region as 
well as far away from wall boundaries. For flow 
regions far away from solid boundaries, the high 
Reynolds number form of the model can be used; 
however, wall-damping functions must be used near 
wall boundaries. 

Speziale (ref. 5) developed a nonlinear version of 
the k-e model that broadened the range of linear 
model application but maintained most of the pop- 
ular features (such as reduction to mixing layer the- 
ories for thin shear flows and the ease of application 


in existing Navier-Stokes solvers without a substan- 
tial increase in the computational time). Speziale 
developed the new model by making an asymptotic 
expansion subject to constraints of dimensional and 
tensorial invariance, realizability, and material frame 
indifference. The model thus obtained was shown 
to yield substantially improved predictions in incom- 
pressible turbulent channel flows and to yield normal 
Reynolds stress differences that give rise to secondary 
flows in square ducts. In the present research work, 
the nonlinear model developed by Speziale was mod- 
ified to include additional terms that contribute sub- 
stantially to the magnitude of Reynolds stresses near 
the wall boundaries. 

For transonic and supersonic flow propulsion ap- 
plications, the local density variation in standard in- 
compressible models does not adequately duplicate 
the experimentally observed reduction in growth rate 
of the mixing layer with increasing convective Mach 
number. However, substantial progress has been 
made in the development of appropriate compress- 
ibility corrections to the transport equation turbu- 
lence models. (See refs. 6 and 7.) These corrections 
resulted from direct numerical simulation of homoge- 
neous compressible turbulence. Notably, Sarkar et al. 
(ref. 6) recognized the importance of including com- 
pressible dissipation in the two-equation turbulence 
model when computing high-speed flows. A sim- 
ple correction was proposed for compressible dissipa- 
tion that can be included easily in the existing two- 
equation turbulence models. The standard model is 
recovered when the model constants for these correc- 
tions are assumed to be zero. 

The objective of this study was to system- 
atically investigate the effect of grid resolution, 
near-wall damping, and various turbulence models 
on the computed flow field. A general-purpose, 
three-dimensional, multiblock Navier-Stokes code 
(described and applied in refs. 8-10) was used in the 
present study. The flow solver contains the Baldwin- 
Lomax turbulence model, a two-equation k-e turbu- 
lence model with various near-wall damping func- 
tions, and a nonlinear stress model for resolution of 
flow-field anisotropies. In addition, the code has a 
built-in performance module to compute quantities 
such as lift, drag, thrust, and discharge coefficients. 
During a typical numerical simulation, these quan- 
tities are constantly monitored to assess the perfor- 
mance of the propulsion system. 

The computed results were compared with pub- 
lished experimental data for flow fields of increas- 
ing complexity. The geometries considered were a 
flat plate (ref. 11), 16° and 24° compression corners 
(ref. 12), a two-dimensional airfoil section (ref. 13), 


and supersonic flow through a square duct (ref. 14). 
The flat plate was selected because it is the simplest 
of all the geometries for which the effects of vari- 
ous near- wall grid spacing, turbulence models, and 
damping functions can be tested. The compression 
corner and airfoil geometries represent the next level 
of flow complexity because these cases contain a sep- 
arated flow region that interacts with a shock.. In__the 
square duct, a secondary flow structure develops per- 
pendicular to the main flow and is mainly attributed 
to the flow-field anisotropy which is not simulated 
by linear (isotropic) models. Depending on the na- 
ture of the flow, either the space- or time-marching 
options in the PAB3D code can be used. Space- 
marching solutions were obtained for the flat plate 
and square duct geometries. Time-dependent options 
were used to investigate separated and transonic flow 
fields (compression ramp and airfoil geometries). 

Symbols 

C cp , C wk constants in Baldwin-Lomax 
turbulence model 


F(n) function in Baldwin-Lomax turbu- 

lence model 

F max maximum of F(n) 

F wake wake function for Baldwin-Lomax 

turbulence model 

fp damping function for k-e equations 

H Heaviside function 

J Jacobian of coordinate 

Q 

transformation, m 

k turbulent kinetic energy, Pa 

L/j. near-wall term for k equation 

L e near- wall term for e equation 

l mixing length for turbulent 

viscosity, m 

M Mach number 

M t turbulent Mach number 

Mto cutoff turbulent Mach number for 

Wilcox model 


a 

duct half-width, m 

Cd,Ce 

model constants for nonlinear 
model 

Cp 

average skin-friction coefficient, 
Jo Cf d x 

Cf 

rr» • ^“W 

local skin-friction coefficient, — 

<?oo 

c p 

pressure coefficient 

Cp 

turbulence viscosity coefficient for 
k-e model, 0.09 

c 

airfoil chord length, m 

Cd 

section drag coefficient 

ci 

section lift coefficient 

D 

duct width, cm 

E, F, G 

convective terms in x, y, and z 
directions 

E.F.G 

flux terms in £, 77 , and £ 
directions 

G„ 

diffusion terms in x, y, and 2 
directions 

e 

total energy, Pa 

Eftic 

total vector skin-friction force, N 

•fkleb 

Klebanoff intermittency factor 


n 

N (n 1 ,n 2 ,n 3 ) 


^max 

P 

Pr 

Prt 

P 

Q 

Q 

Q 

R c 

Rt 

Rx 

S 

S k 


normal distance from wall, cm 
unit normal vector 


law-of-the-wall coordinate, 

TiyJ pwi~w np w u T 

fJ'W Pw 

normal distance from wall at F m3X 
location, cm 

production term for k-e equations 
Prandtl number, 0.75 
turbulent Prandtl number, 0.9 
pressure, Pa 
conservative variable 


= Q 

J 

heat flux 


Un 

cell Reynolds number, R c = 


turbulent Reynolds number 

Reynolds number based on stream- 
wise distance from plate leading 

j UocPocX 

edge, 

Moc 

source term for Navier-Stokes 
equation 


source term for k equation 
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Se 

source term for e equation 

T 

temperature, K 

t 

time, sec 

u,v,w 

velocity components in f , 77 , and £ 
directions, m/sec 

u , v , w 

velocity components in x, y, and z 
directions, m/sec 

u + 

law-of-the wall coordinate, 
u _R C 
u T n+ 

U T 

friction velocity, </— 

V P 

x,y,z 

spatial coordinates, cm 

T 

compressibility correction factor 

7 

ratio of specific heat, 1.4 

AA 

incremental cell face area, m 2 

6 

Kronecker delta and boundary 
layer thickness 

€ 

turbulent energy dissipation 

e 

boundary layer momentum 
thickness, cm 

K 

Von Karman constant 

P 

dynamic viscosity 
coefficient, m 2 /sec 


generalized coordinates as func- 
tions of x, y, and t 

p 

density, kg/m 3 

a k 

diffusion coefficient for k equation 

&€ 

diffusion coefficient for e equation 

T ij 

Reynolds stress components, where 
i and j represent x, y ) or z 

r w 

wall shear stress, 

an \w 

oj 

vorticity 

Subscripts: 


cros 

matching point for inner and outer 
boundary layer regions 

e 

edge of boundary layer 

i 

inner 

max 

maximum 


mm 

minimum 

0 

outer 

w 

wall 

x, y,z 

x, y, and 2 derivatives 

00 

free stream 

Superscripts: 


L 

laminar 

T 

turbulent 

Abbreviations: 


k-e 1 

k-e model, Jones and Launder wall 
damping 

k-e2 

k-e model, Van Driest wall 
damping 

k-e 3 

k-e model, Speziale et al. wall 
damping 

RMS 

residual mean square 

WF 

wall function 

Experimental Configurations 


Flat Plate 


The supersonic data for this model were ob- 
tained for an insulated flat plate tested in the NASA 
Ames Research Center 6 -inch Heat- Transfer Tun- 
nel. The model was a flat surface 40.64 cm (16 in.) 
long and had a lower surface leading-edge chamfer 
of 15°; it spanned the width of the tunnel. (See 
ref. 11 .) The leading edge was rounded to a radius 
of 0.0076 cm (0.003 in.). The static pressure ori- 
flees were 0.0343 cm (0.0135 in.) in diameter and 
were placed 2.54 cm (1 in.) from the side edges of 
the plate. Transition from laminar to turbulent flow 
was forced near the leading edge by a strip of lamp- 
black 1.27 cm (0.5 in.) wide placed at the leading 
edge of the plate. Data were obtained at a free- 
stream Mach number of 2.5 and a tunnel total pres- 
sure of 204 kPa (30 psia). The Reynolds number 
based on the distance measured from the flat-plate 
leading edge ranged from 2.1 x 10 6 to 6.2 x 10 6 . 

Boundary layer measurements were made over a 
survey area 10.16 cm (4 in.) to 20.32 cm (8 in.) 
from the plate leading edge. The boundary layer 
measurements were performed with a total-pressure 
probe. The centerline of the probe was 0.0165 cm 
(0.0065 in.) above the surface when the probe was in 
contact with the plate surface. The probe had a rect- 
angular external dimension of 0.2032 cm (0.080 in.) 
by 0.0330 cm (0.013 in.). 


3 


The free-stream Mach number varied by no more 
that 2 percent in the test section of the tunnel. That 
variation resulted in a skin-friction coefficient error of 
less than 2 percent; the total instrumentation error 
for the skin-friction coefficient was ±2.5 percent. 

Compression Ramps 

The compression ramp investigations were con- 
ducted in the 8 x 8-in. Supersonic Blowdown Tun- 
nel at Princeton University. (See ref. 12.) The tun- 
nel can provide a test duration of 30 sec to several 
minutes at stagnation pressures of 5.1 to 50.7 MPa 
(50 to 500 atm) over a Mach number range of 2.84 
to 2.95, depending on which test sections are uti- 
lized. All of the models were tested in nearly adia- 
batic wall conditions. The particular data used for 
this study were obtained at a Mach number of 2.85 
and a stagnation temperature of 262 K (472 R). The 
corresponding free-stream unit Reynolds number was 
0.64 x 10 6 cm' 1 (1.6 x 10 6 in" 1 ). The solid brass 
model consisted of a short upstream flat segment 
joined to a ramp that was 15.24 cm (6 in.) long by 
16° or 24° compression corners. Sidewall fences were 
attached to each side of the ramp model to lessen the 
influence of the tunnel wall boundary layer on the 
compression corner flow. The boundary layer probe 
used for this study had a 0.0178-cm (0.007-in.) flat 
tip with 0.0076-cm (0.003-in.) orifices. 

Pitot pressure tubes were used to make detailed 
flow- field surveys upstream of the compression ramps 
to determine the approaching flow properties. The 
displacement and momentum boundary layer thick- 
nesses determined from these measurements were 
subsequently used to determine the computational 
inflow boundary conditions. The estimated errors 
were ±5 percent in the streamwise velocity compo- 
nent and ±10 percent in the skin- friction coefficient 
determined from the Preston tube measurements. 

Subsonic Airfoil Section 

The airfoil test case has an RAE 2822 airfoil con- 
tour (ref. 13) and is a subcritical design section with a 
trailing edge thickness of zero. The airfoil is 12.1 per- 
cent thick and is designed for a lift coefficient of 0.56 
at a free-stream Mach number of 0.66. The data were 
obtained for a model with a 0.61 m (2 ft) chord and 
that spanned 1.83 m (6 ft) during an experiment con- 
ducted in the RAE 8 x 6- ft Transonic Wind Tunnel. 
(See ref. 13.) The tunnel is a continuous, closed cir- 
cuit type that operates at a stagnation pressure range 
of 10 to 355 kPa (1.5 to 50 psi) with an average stag- 
nation temperature of 307 K (552°R). Surface static 
pressure data, wake pitot and static pressure data, 


and boundary layer pitot and static pressure data 
were obtained for a variety of conditions with the 
transition fixed. Additionally, the oil-flow visualiza- 
tion technique was used to observe flow separation. 
An extensive description of the tunnel flow condition, 
wall interference, and instrumentation can be found 
in reference 13. 

Square Duct 

The experiment was set up in a continuous flow, 
open-circuit wind tunnel with a test section in the 
form of a square duct 50.8 cm (20 in.) long made of 
Plexiglas 1 material; it had a cross section of 5.08 cm 
(2 in.) at the duct inlet. A square brass constant- 
area duct with 2.54 cm (1 in.) on each side was 
placed within the outer Plexiglas duct. This “duct 
within a duct 1 ’ configuration was built to ensure a 
clean starting condition for the inner duct flow by 
allowing the distorted flow that develops along the 
side walls of the outer duct nozzle to be bypassed 
through the annular space between the inner and 
outer ducts. The experiments were conducted at a 
free-stream Mach number of 3.9. The total pressure 
and total temperature at this location were 276 kPa 
(40 psi) and 300 K (540°R), respectively. A circu- 
lar pitot tube with an outside diameter of 3.05 mm 
(0.12 in.) was used to obtain total pressure profiles at 
three streamwise locations; a flattened-tip probe with 
outside dimensions of approximately 2.5 x 6.6 mm 
(0.1 x 0.26 in.) was used to obtain boundary layer 
measurements. Wall shear stresses were calculated 
from measurements obtained with several differently 
sized Preston tubes resting on the wall. The experi- 
mental uncertainty in measuring the skin-friction co- 
efficient is about 10 percent. Other pertinent details 
of the experimental setup and instrumentation are 
given in reference 14. 

Theoretical Formulation 

The governing equations of the Reynolds- 
averaged Navier-Stokes formulation include the con- 
servation equations for mass, momentum, energy, 
and the equation of state. In the present study, the 
perfect gas law is chosen to represent the properties of 
air. For flow that contains turbulence, the Reynolds 
stresses are modeled using the eddy viscosity concept. 

Turbulence models are essential for the realis- 
tic simulation of aerodynamics in the high Reynolds 
number regime. Two turbulence models were used 


1 Plexiglas: Rohm and Haas Co., Inc., Philadelphia, 

Pennsylvania. 
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for the current study: the Baldwin-Lomax algebraic 
turbulence model (ref. 15) and a two-equation k-e 
turbulence model that follows the formulation of 
Jones and Launder. (See ref. 1.) The Navi er- Stokes 
equations and the mathematical representations of 
the relevant turbulence models are described in the 
following sections. 


Navier-Stokes Equations 


The mass, momentum, and energy conservation 
equations of the Reynolds-averaged Navier-Stokes 
equations can be written in terms of generalized 
coordinates and in a conservative form as follows: 


dQ dE dF dG 
dt + d£ + d C + 9^ -S 


( 1 ) 


where S = {0, 0, 0, 0, 0} T for laminar or algebraic 
turbulence modeling. However, if a turbulence ki- 
netic equation is used, the source term for the en- 
ergy equation is replaced by the source term -S k . 
In equation (1), 


e = 1[^e + ^f + 6G] 

-K*E + fcF + &G] 

F = ±[(CsE + Cs,F-K*G) 

— (CrE„ + £yF v + 

G = j[(Vx E + r) y F + rj z G) 

— (T]xF V + TJyFy -I- 77 2 Gi,)] 


f p 

pu 

Q pv 

pw 


e 


{ pw t 

pvw > 
pw 2 + p 
(e 4- p)w J 


E v ={ 


‘zy 

T xz 

-Qx + UT XX + VT xy -(- WT X 
0 


Ft; = < 




'xy 

T yy 

T yz 

-q y + UT X y + VTyy + WTy z 

0 

T xz 

Tyz 

T Z Z 

-q Z + UT XZ + VTyz 4- WT ZZ J 


-qx = 


~Qy — 


f* 


fi T \ da 2 

7 — 1 ^ Pr ' Prt J dx 
1 U L 

7 — 1 ( Pr Prt J dy 


+ 


Qz 




7—1 \Pr 



In these equations, p is the density; u, v, and w are 
the velocity components in the x, y, and z directions, 
respectively; e is the total energy per unit volume; 
the pressure p is related to e by 


P = ( 7 ~ 1) e — + v 2 + w 2 ^ 


( 2 ) 


and, for example, r xy = t^ + t^. The laminar shear 
stress r^y may be expressed in the following forms: 


r r i 

\ pu*+p 

E = < puv ► 
puw 

l (e + p)u > 


F = < 


pv 

puv 

pv 2 +p 


> 


pvw 

l (e +p)v J 




dv 

dy 


dw 

dz 


’xy 



( 


du <9tA 
dy ^ dx ) 


yy 




dw 

dz 


du 

dx 


r yz — f 1 ( 


dv 

\dz 
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r 2 t ( dui du dv\ 
t " = 3 M Vlh ~ dy) 



The forms of the turbulent shear stress rj y for various 
turbulence models are described in the sections that 
follow. In the PAB3D code, all streamwise deriva- 
tives of the Reynolds stress (d/d£) have been omit- 
ted from equation (1) for computational economy. 
This formulation is known as the simplified Navier- 
Stokes equation. The remaining cross-stream deriva- 
tives can be numerically implemented either coupled 
or uncoupled. The uncoupled derivative is the thin- 
layer Navier-Stokes approximation. In all the results 
in this paper, the derivatives are uncoupled. 

Baldwin-Lomax Turbulence Model 

The Baldwin-Lomax model (ref. 15) is an alge- 
braic two-layer turbulence model in which the tur- 
bulent eddy viscosity is evaluated as follows: 


For the inner layer, 


(m t )» = P i 2 M 

where the mixing length for turbulent viscosity 
l = 0.4njl — exp(— n + /A + jj 
and the vorticity 


( 3 ) 


( 4 ) 


w = y/(uy - U X ) 2 + (v 2 - Wy ) 2 +(w x - u z ) 2 (5) 
For the outer layer, 

(/I^ 1 )^ = O.OlGSC^cpP^wake^klebC 7 ^) (®) 

where F w ^ke smaller of ri max F max or 

C’ w i c 77, max {'u+v 2 +w 2 ) max F } max ■ The term n max is the 
value of n that corresponds to the maximum value of 
the model function F, F ma xi where 


fjiF — (t^)i { n — n cros) 

f? = (jl T )o ( n n cros) 

where n is the normal distance from the wall and 
n C ros is the smallest value of n at which magnitudes 
of the viscosities at the inner i and outer o boundaries 
and {(jl T ) 0 are equal. The turbulent stress is 
determined from 

( du z duj \ 2duk 
\dxj + dx t ) 3dx k lJ 



F(n) = n|u>| [l - exp(-n + /"4 + )] (7) 

and the Klebanoff intermittency factor F kleb is cal- 
culated by 


Fkleb 


i+5.5f 

V TT-max / , 


-1 


( 8 ) 


The values of the constants appearing in equa- 
tions (3)— (8) are listed in reference 15 as A + = 26, 
Cep — l-6j C wk = 0.25, and C k i eb = 0.3. 


Two-Equation k~E Turbulence Model 


The Jones and Launder (ref. 1) formulation for the two-equation turbulence model uses k and e as the 
principal variables. A modified form of the original Jones and Launder model is used in this study. This 
modified formulation is fully three dimensional, and the governing equations are written in a conservative form 
in terms of generalized coordinates. The governing equations can be cast in the same form as the Navier-Stokes 
equations (eq. (1)) with the following new definitions of the dependent variable and source terms: 
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S £ — Ci P— — C 2 p— + L e S k — P p{ 1 + r)e + Lf. 



p_ T du T dv T dw 

xx dx yya y zz dz 


+ T, 


T ( du dv ^ 


xy V dy dx 


M , T f \ r / dw du \ 

1) y *\dz + d y y Tzx Xte + ^) 


or is expanded to 


where 


p _ n T(( du , dv\ 2 fdv dw\ 2 (dw du\ 2 

F ~' 1 (l % + &J + (& + W + (fe + aj) 

4 - o\{ —Y -L ( 9v \ 2 > f dw \ 2 ] 2 (du dv dw\ 2 ) 

LV^xJ \dy) + \dz) J 3 \dx + dy* ) 


2 ( du dv dw\ 

j'*!.* + 5 + & j 




1. 2 

<>7 




Mfc = ^ C M = 0.09/^ 


Cj — 1.44 C 2 = 1.92^1 — 0.3 exp 
CT e = 1.3 cr k = 1.0 


MIC 


t r 

T zj = M 


<9tij 

3x, + <9 x 7 


2 3m. 


3 dx; 6ij 


2 pk^ij 


( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


where * and j represent x, y, or z. The compressibility correction T and the damping function f u are described 
next. M 


Compressibility correction functions for k-e 
model. High-speed turbulent flows have different 
characteristics than low-speed flows (incompressible 


flows). For example, the rate of spread of the shear 
layer in high-speed flows is much slower compared 
with that of low-speed flows. Several corrections for 
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this effect have been developed in the last few years. 
The two most widely used compressibility corrections 
are given below. 

Sarkar et al model (ref. 6): 

r = M t 2 (i6) 

Wilcox model (ref. 7): 

r = {mJ — MfQ^H(M t — M t o) (17) 

where H (x) is the Heaviside step function; Mt is the 
local turbulent Mach number defined as \fkf a, in 
which a is the local speed of sound; and T = 0 for no 
compressibility correction. 

Damping and wall functions for k-e model. Solv- 
ing wall-bounded flows requires the use of damping 
or wall functions. These functions adjust the turbu- 
lence viscosity near solid surfaces. The details of the 
damping and wall functions are given as follows. 

Damping functions . The damping function 
adjusts the turbulent viscosity through the term C^. 
Far from the wall, fp = 1; at the wall, = 0. 
Three damping function forms were investigated in 
this study. They are referred to as ft-el, k-e 2, and 
k-e 3 and are defined as follows. 

The k-e 1 (Jones and Launder) form (ref. 1): 
ffi = exp 

The k-e 2 (Van Driest) form (as used by Nagano 
and Hishida, ref. 16): 


3.41 


1 + (ft/60) J 


(18) 


grid spacing normal to a wall with n + > 5. One way 
to achieve greater accuracy for grids with n + < 50 
is through the use of the simple wall function (WF), 
although the use of this function is limited to the 
calculation of attached flows. The concept estimates 
the wall shear stress from law-of-the-wall coordinates 
n + and of the first cell from the wall (he., a sur- 
face cell). The equivalent turbulence viscosity at the 
wall is subsequently calculated from the estimated 
shear stress. (In the standard two-equation turbu- 
lence model, the turbulence viscosity at the wall is 
normally set to zero.) The wall function approach 
can be used to relax the restriction on n + , which 
permits use of values up to 50. This approach typ- 
ically speeds solution convergence rates by reducing 
the total grid count to describe the problem. The 
following steps describe the procedure for the wall 
function approach. 

First, evaluate the surface cell Reynolds number 

33 D Pwun 

H c — 

Pw 

The Reynolds number R c can also be written as a 
function of the two nondimensional parameters u + 
and as R c = Thus, 


^ (21) 

A normalized velocity profile that relates and 

near a solid wall is described in many fluid 

mechanics textbooks. In a turbulent flow along a 
wall, Von K&rman and others have suggested that 
the flow should be divided into three zones governed 
by the value of n + . Therefore, 

u + = /(n + ) (22) 


U = i - ex P -jr ( 19 ) 

Also, the quantity 2 in the near- wall term L e (eq. (9)) 
is replaced with 1 — f )L - 

The fc-e 3 (Speziale, Abid, and Anderson) form 
(ref. 17): 



On any solid surface, the dissipation is set equal to 
Lfc. Then, L e and L* are set to zero. 

Simple wall function form. Neither the Baldwin- 
Lomax nor the two-equation turbulence model (even 
with damping functions) is capable of producing ac- 
curate aerodynamic predictions based on minimal 


In turbulent flow along a wall, the flow may 
be divided into three regions. First is the laminar 
sublayer in which the viscous stress is much greater 
than the Reynolds stress: 

ti + = n + (n + < 5) (23) 

The region immediately above the laminar sublayer 
is called the buffer zone. In this zone, the Reynolds 
and viscous stresses are of the same order: 

u + = 11.5 log 10 n + - 3.05 (5 < n + < 30) (24) 

Farther from the wall, in the turbulent zone, the 
Reynolds stress is much greater than the viscous 
stress. Thus, 

u + — 5.75 log, o n~^ + 5.50 (30 < n^~ ) (25) 
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In the solution to equations (21) and (22), n + is 
evaluated and the value of t w can be specified at the 
wall as follows: 


T w — 


[(fan + )/n] 2 

Pw 


(26) 


Then, because t w = (f* L + H T )$%\ W , a corrected tur- 
bulent viscosity /ij at the wall can be obtained. This 
corrected viscosity replaces the turbulent viscosity in 
all the transport equations. 


anisotropies. The nonlinear terms added to the stan- 
dard model are important for the prediction of sec- 
ondary flow in a square duct. The secondary flow de- 
velops in the cross-flow planes due to the cross-stream 
gradients of the Reynolds stress. For linear mod- 
els the cross-stream gradient of the Reynolds stress 
difference is small; therefore, anisotropic flow fea- 
tures do not develop. However, the nonlinear mod- 
els produce gradients of sufficient magnitude to de- 
velop the anisotropic or secondary flow feature. For 
Cd = Ce = 0, the linear model is recovered. 


Nonlinear k-e model. The nonlinear k-e model 
is obtained by adding quadratic terms to the linear 
model rf- to treat the mean velocity gradients. In 
equation (27), the first two terms on the right-hand 
side correspond to the linear model and the addi- 
tional terms represent the contribution due to non- 
linear effects. These changes to rj- in E v , F^, and 

as well as P are as follows: 


tT.= 


-it 2 / T d ujc \ 


+ 2 ,?D t j+C D Cl^ 


x ( DimDmj ^DmnDmn&ij 


+ C E Cl~- (Z)y - |jD mn <5y) 


2 L k(dVk \ 2 

+ ^ ^ (n,m) (27) 


W v (n,m) - -Sij - 6 in S jn + 46 im 6. 


u jm 


(28) 




- duj — 

dx k J 


duj 

dx k ‘ 


where 


if duj du ? 


Dij — - 


l i ^“7 

2 \ dxj dxj 


Cj) and Ce axe model constants, and Cd = 
Ce — 1.69. 


This new nonlinear model differs from the origi- 
nal model developed by Speziale (ref. 5) by the addi- 
tion of the last term on the right-hand side of equa- 
tion (25). (See development in ref. 18.) Because both 
k and e vary rapidly near solid boundaries, the addi- 
tional term contributes significantly to the near-wall 


Performance Method 

The performance method (refs. 19 and 20) obtains 
body forces through the application of the momen- 
tum theorem to a control volume that surrounds the 
model. The choice of surfaces over which the integra- 
tion of forces is performed provides several options for 
calculation of the momentum and pressure forces on 
the model. The method used for this investigation 
integrates the mass flux and pressure forces over the 
model with 

F = £>U(U • N) + (p - Poo)N] A A + F fric (30) 

where F is the total vector body force, AA is the 
area attributed to the cell face, and N is the unit 
normal vector of the cell face. The static pressure 
force on a solid wall is calculated by extrapolating the 
cell-centered static pressure to the wall surface and 
by assuming a zero velocity at the wall. The term 
U • N for solid walls vanishes as solution convergence 
is obtained. 

Skin friction is calculated for the solid wall bound- 
aries of the control volume. The viscous stress ten- 
sor used to determine the skin-friction force is cal- 
culated with only the velocity derivatives normal to 
the surface. The velocity gradients are determined 
by a two-point difference. The first velocity is a zero- 
magnitude vector positioned on the surface. The sec- 
ond velocity is the velocity at the cell center. The 
local shear stress tensor is constructed from the nor- 
mal velocity gradients multiplied by the local viscos- 
ity. The viscosity was determined from Sutherland’s 
formula (ref. 21) and used the static temperature at 
the local cell center. 

Method of Solution 

The simplified Reynolds-averaged Navier-Stokes 
equations and the associated turbulence models have 
been implemented in the computer code PAB3D. As 
mentioned previously, the numerical code has the op- 
tion for either space- or time-marching solutions. In 
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particular, the space-marching option is well suited 
for supersonic flows with no embedded subsonic re- 
gion or for flows in which pressure gradients are 
practically absent. The space-marching algorithm in 
the PAB3D code is also much more robust than the 
parabolized marching procedure. For cases in which 
the space-marching scheme criterion is met, the solu- 
tion is as accurate as those obtained with the time- 
dependent algorithm. However, the space-marching 
procedure requires significantly less computer time 
(approximately a factor of 20 less) than does the 
time-dependent procedure. 

Computational Procedure for Navier-Stokes 

Equations 

The solver for the Navier-Stokes equations has 
been implemented in the PAB3D code with three dif- 
ferent numerical schemes: the Van Leer flux- vector 
splitting scheme, the Roe flux-difference splitting 
scheme, and the space-marching scheme which is a 
modified version of the basic Roe scheme. These 
schemes are basically implicit and upwind and are 
constructed by using the finite volume approach. 
Only the inviscid portion of the flux vectors E, F, 
and G are subjected to the splitting and upwind pro- 
cedures. The diffusion terms of the Navier-Stokes 
equations are centrally differenced. A detailed de- 
scription of the mathematical formulation for these 
schemes can be found in reference 8. 

The flux-vector or flux-difference splitting is used 
in all three computational directions. The updated 
solution at each iteration is obtained by using an 
implicit procedure in the rj and £ mesh planes at 
each constant £ value. The relaxation procedure in 
the £ direction consists of a forward and backward 
sweep. This particular implementation strategy has 
an important advantage: because the metrics for 

the implicit procedure are only required for up to 
three planes, they are not stored for the entire grid 
domain. Instead, they are recomputed one plane at 
a time at the advancing front of the prevalent sweep 
direction. This approach requires less memory for the 
intermediate data structure. Typically, 20 words of 
memory are required for each grid point for moderate 
to large mesh sizes. 

For a general time-dependent solution that uses 
the Van Leer or the Roe scheme, each iteration 
count contains a forward and backward sweep in the 
i direction, with one step of an implicit update of the 
solution in each of the cross planes. For several super- 
sonic and subsonic flow conditions, the numerical 
scheme of Roe can be further simplified into a space- 
marching method as follows. The inviscid terms 


in the Navier-Stokes equation are discretized as an 
approximate Riemann problem. The interface flux 
in the streamwise direction is determined by separate 
terms that depend on the quantities on the upstream 
and the downstream sides of the interface. For fully 
supersonic or subsonic flow with a small pressure 
gradient, the information can travel only in the flow 
direction and is carried by the terms on the left- 
hand side. For these flow problems, the upstream 
effect carried by the terms on the right-hand side 
can be ignored when compared with the streamwise 
influences. A solution is obtained by performing 
sufficient implicit iterations in each plane until the 
convergence criteria are met, A solution in the 
entire computational domain is established in a single 
forward sweep. All solutions in this paper were 
obtained with either the standard or space-marching 
version of the Roe scheme. 

Computational Procedure for k-£ Equations 

The governing equations of the two-equation tur- 
bulence model are written as a pair of coupled trans- 
port equations in conservative form. In principle this 
model could be implemented with the Navier-Stokes 
equations as a set of seven coupled equations, or the 
model could be a separate implementation that is un- 
coupled from the Navier-Stokes equations. The fully 
coupled approach would result in an increase in the 
computational requirements and numerical stiffness. 
For this study, the k-e equations are implemented un- 
coupled from the Navier-Stokes equations and from 
each other. Although the stiffness remains, it is al- 
leviated to some degree by solving these two equa- 
tions with a much smaller Courant-Friedrichs-Levy 
(CFL) number (usually 0.25 of the CFL number of 
the Navier-Stokes equations). The potential differ- 
ences in the development of the flow over time and 
turbulence variable sets have not noticeably affected 
the convergence rate or the quality of the solutions. 

The governing equations for the nonlinear model 
are the same as for the linear model except for 
the differences in the expressions for the Reynolds 
stresses. Because the additional nonlinear terms in 
the Reynolds stresses are not large, they are treated 
simply as added source terms in the code. However, 
the variables ti, u, it;, fc, £, and their first deriva- 
tives are already calculated for the linear model. The 
new variables that need to be calculated are the sec- 
ond derivatives. When the first and second deriva- 
tives are known, the nonlinear contribution of each 
component of the Reynolds stresses can be obtained. 
The nonlinear model required only 2 percent addi- 
tional computational time at each time step and con- 
verged somewhat more slowly than the linear model. 
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This result is not surprising because the nonlinear 
k-e model contains Reynolds stress relaxation terms 
that are dispersive rather than dissipative. 

Multiblock Structure and Boundary Conditions 

The PAB3D code uses a general multiblock grid 
topology to handle complex configurations. One-to- 
one, multiple- to-one, or general patched interfaces 
between the blocks are accepted by the code. An 
important feature of the PAB3D code is the partition 
concept in the streamwise direction. If different 
viscous stress models are employed within a block, 
the length of the block can be partitioned by means 
of the starting i index for each viscous stress model. 

The boundary conditions often include inflow, 
outflow, free stream, solid walls, and geometrical 
symmetry. Five types of inflow and outflow boundary 
conditions are provided: Riemann characteristics, 
fixed inflow total temperature and pressure, com- 
pletely fixed inflow parameters, constant pressure for 
subsonic outflow, and extrapolation for supersonic 
outflow. The Riemann characteristics boundary con- 
dition is used at free-stream boundaries. On a solid 
boundary, either a no-slip or an inviscid-slip bound- 
ary condition can be specified. Finally, the sym- 
metry boundary conditions include mirror imaging 
across a plane and polar symmetry around an axis 
in any direction. A universal high-order symmetry 
boundary condition is used that was developed by 
Abdol-Hamid and Pao. (See appendix.) A logically 
simple control structure is required to direct the code 
execution for dimensions of the zones and blocks, 
solver options, connections between blocks, bound- 
ary conditions, time-stepping requirements, and tur- 
bulence models. 

For the turbulence transport equations, either 
zero-order extrapolation or free-stream values are 
used for k and e along the outer boundaries. If 
the flow is outgoing along the outer boundary, zero- 
order extrapolation is used. If flow entrainment 
is involved, then free-stream values are used along 
the outer boundaries and the values of k and e 
are set such that a preselected nominal turbulence 
intensity is achieved. For all the results in this 
paper, this level is 1 percent of the flow velocity. For 
the inflow condition, laminar or algebraic turbulence 
model solutions were computed for the first few 
planes of the upstream blocks. In that flow region, 
the inflow profile for k takes the same shape as 
the vorticity profile except that it is multiplied by 
a specified value of maximum turbulence intensity. 
(See fig. 1.) When the k profile is known, the e profile 
is obtained based on the hypothesis that production 


equals dissipation. These values of k and e are used 
as boundary conditions for the k-e blocks. On any 
solid surfaces in which the k-e 3 model is applied, the 
dissipation e w is set equal to the value of L k (eq. (10)) 
and k is set to zero. For k-e2 and k-e2 , e w = k w = 0. 

Results and Discussion 

We selected the well-established properties of flow 
over a flat plate to calibrate the different forms of the 
two-equation k-e model utilized in this report. Mach 
numbers of 2 and 0.4 were selected for comparisons. 
Figure 2 shows law-of- the- wall solutions (u + = u* /u r 
versus y+) that use the k-e 1 model for subsonic and 
supersonic flows at R e = 30 000, where u* is defined 



Similar results were produced for both subsonic and 
supersonic flows. Then, the supersonic flow case was 
computed based on the Baldwin-Lomax turbulence 
model. (See fig. 3.) Both turbulence models agreed 
well with the theoretical curves. The k-e 1, k-e2 , and 
k-e 3 turbulence models are compared (fig. 4) and all 
three models produce very similar results; that is, 
all three models predicted the Von Karman constant 
of k — 0.41 within less than 2 percent. As expected, 
the different forms of the turbulence models give very 
similar predictions for this attached flow case. 

Supersonic Flow on Flat Plate 

In many cases, skin-friction drag represents a sig- 
nificant portion of the total drag of a supersonic ve- 
hicle. Accurate prediction of skin friction is essential 
for CFD applications in design and analysis. An in- 
sulated flat plate (ref. 11) operated at a Mach num- 
ber of 2.5 over a Reynolds number range of 2.1 x 10 6 
to 6.2 x 10 6 is modeled in the present study. Pre- 
dicted velocity and average skin-friction predictions 
are compared with the experimental data. Both the 
Baldwin-Lomax and the k-e turbulence models are 
used to predict the aerodynamic characteristics of 
this flow. 

In the present analysis, a grid distribution of 
71 x 2 x 81 was used. The plate was 12 in. long, 
and the Reynolds number at the back of the plate was 
taken to be 6.2 x 10 6 . First, the number of grid points 
was fixed at 81 in the k direction (normal to the 
wall) and the value of n + (nondimensional distance 
normal to the wall of the first grid) varied from 0.5 
to 10. Figure 5 shows the streamwise velocity profile 
(at x — 6 in.) versus y/0 , where 6 is the boundary 
layer momentum thickness. Both turbulence models 
predict reasonably well the overall boundary layer 
velocity development for < 10. However, the 
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k-e 1 model kept a very well-defined boundary layer 
development at n+ = 10. On the other hand, the 
Baldwin-Lomax turbulence model gives a boundary 
layer development for n + = 10 similar to the one 
produced from a laminar solution. Figure 6 shows 
the average skin-friction predictions compared with 
experimental data as a function of Reynolds number. 
The Reynolds number is calculated based on the free- 
stream condition and the distance along the plate 
from the leading edge. The k-e model average skin- 
friction prediction agreed well with experiment for 
n* < 5. However, the Baldwin-Lomax turbulence 
model gave poor predictions of Cp for n + > 2.5. 

We performed a second set of computations by 
fixing = 5 and changing the number of grid points 
in the normal direction from 21 to 81 with the k-s 
turbulence model. Figure 7 shows the comparison 
of experimental data for a Baldwin-Lomax solution 
computed with 81 grid points in the normal direction 
and the k-e turbulence models solutions with varying 
numbers of grid points. The k-e model produces 
similar results for all grid distributions computed. 
A similar result was obtained for the average skin 
friction on the plate. Thus, the distance of the 
first grid point from the wall is a very important 
parameter in predicting boundary layer flows, but 
the number of grid points is not important. 

As shown in figure 6(a), the two-equation k-e 1 
turbulence model fails to predict the wall skin friction 
for n + > 5. Figures 8 and 9 show the improvement 
obtained by using the simple wall function in k-e 1 
to predict the streamwise velocity variation and skin 
friction at = 10 and 50, respectively. 

The flat-plate test case was selected because it 
is the simplest of all the geometries for which the 
effects of various near-wall grid spacing, turbulence 
models, and damping functions can be easily tested. 
The results obtained with the k-e turbulence model 
were much more accurate and consistent than those 
for the Baldwin-Lomax turbulence model for a wider 
range of (height of first grid line normal to 
solid surface). The n ^ value is a very important 
parameter for accurate prediction of boundary layer 
development in attached and separated flows. The 
number of grid points normal to the solid surface is 
not as important as a correct value of n + . 

16° and 24° Compression Ramps 

The compression ramp is a very simple test case 
but contains complex aerodynamic characteristics. 
A flow of Mach 2.85 over 16° and 24° ramps has 
been computed and compared with data by Settles, 


Yas, and Bogdonoff. (See ref. 12.) The free- 
stream unit Reynolds number was 0.64 x 10 6 cm" 
(1.6 x 10 6 in^ 1 ). For the 16° compression comer, 
both the Baldwin-Lomax and k-e solutions were ob- 
tained. Only k-e solutions were obtained for the 
24° compression ramp because a steady-state solu- 
tion could not be obtained from the Baldwin-Lomax 
turbulence model. Figures 10 and 11 show the com- 
puted wall pressures and local skin- friction distribu- 
tions and comparisons with the experimental data 
for ramps at 16°. The calculations were done with a 
two-block configuration. The first block employs the 
Baldwin-Lomax turbulence model with 10x2x101 
grid points and sets the initial boundary condi- 
tions for the values of k and e for the downstream 
second block. The second block spans the region 
from upstream of the separated flow region to the 
outflow boundary. This block had the dimensions 
152 x 2 x 101 and used the k-e 1 code option to simu- 
late the flow field. Mesh sequencing is used to assess 
the effect of grid refinement on the computed results 
for the 24° compression ramp. The grid is setup such 
that the first index refers to the number of points in 
the j direction and the last index refers to the num- 
ber of points in the k direction. The base grid, the 
half grid in the j direction, and the half grid in the 
j and k directions were selected to investigate the 
effects of grid resolution on the k-e solutions. Only 
the fine grid (base grid) solution is shown for the 16° 
ramp. In the simulation of the 16° ramp (fig. 10(a)), 
good agreement was obtained in the prediction of the 
surface pressure distribution when both turbulence 
models were used. However, the skin- friction distri- 
bution produced different results. (See fig. 10(b).) 
For the 24° compression ramp, the computed pres- 
sure distribution was insensitive to the effect of grid 
refinement as shown in figure 11(a). However, the 
skin-friction results showed some sensitivity to the ef- 
fect of grid refinement. (See fig. 11(b).) Good agree- 
ment was obtained for skin friction on the upstream 
portion of the ramp. The experimental data on the 
24° ramp indicate a massive separation. The region 
between x/6 = -1.5 and 1.0 where no experimental 
data were taken indicates the separation. The k-e 1 
model was able to predict the attachment point accu- 
rately at x/6 « 0.03, but it predicted the separation 
point slightly downstream of the experimental data 
at x/6 ~ —0.15. 

Airfoil Configuration RAE 2822 

Experimental data for the RAE 2822 (case 10 
at R = 6.2 x 10 6 ) two-dimensional airfoil section are 
available in reference 13. The experimental data 
used for this comparison were at a free-stream Mach 
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number of 0.75 and at an angle of attack of 3.19 
in tunnel coordinates. The computational results 
were performed at 2.8° to compensate for wall in- 
terference; the computation also included an angle- 
of-attack correction equation from reference 13. Six 
blocks are used in the present investigation and to- 
taled 31 000 grid points. Two of these blocks were at 
the leading edge at x/c < 0.03, where the flow was 
assumed to be laminar as in the experiments. The 
outer boundary extent was 8c, and no point vortex 
correction was included in the far field. The base 
grid, the half grid in the k direction, and the half grid 
in the k and j directions are selected to investigate 
the grid effect on the k-e solutions. Figure 12 shows 
the comparisons among the three different grids and 
the experimental data. No significant differences are 
apparent between the half grid in the k direction and 
the base grid in the predicted pressure, velocity, and 
skin friction. 

In the rest of this section we use the base grid 
for comparisons with the experimental data. The 
value for the first grid height in the computational 
mesh was n+ < 2. Figures 13(a)-13(c) show the 
Baldwin-Lomax and k-e 1 turbulence model predic- 
tions compared with the experimental data. In fig- 
ure 13(a), both models predict a similar pressure co- 
efficient C p distribution, although k-e 1 predicts the 
shock location slightly farther forward, which is in 
better agreement with the experiment. Figure 13(b) 
presents the local skin-friction predictions compared 
with the experimental data. The skin friction is nor- 
malized by the local boundary layer edge condition. 
Both models gave similar predictions for skin fric- 
tion upstream of the shock. However, the Baldwin- 
Lomax turbulence model gave a much better pre- 
diction of the skin-friction data point downstream 
of the shock location. In general both models over- 
predicted the local skin friction to x/c < 0.6. Fig- 
ure 13(c) shows the experimental stream wise velocity 
profile at x/c = 0.9 compared with predictions from 
both turbulence models. The k-e 1 model provided a 
more accurate velocity profile than did the Baldwin- 
Lomax prediction. Different wall-damping function 
forms affect C p , the velocity profile, and the local 
skin-friction prediction as shown in figure 14, which 
compares fc-el, k-e 2, and k-e 3. If n+ < 25 is re- 
laxed (fig. 15), again the k-e 1 model (or any of the 
other forms) significantly overpredicts the local skin- 
friction values. However, use of the simple wall func- 
tion brings the local skin friction within the range 
of the other models (predicted at n + < 2). All wall 
pressure distribution predictions with different wall- 
damping functions in the k-e model and the simple 
wall function are summarized in figure 16. The k-e 1 


model (n + = 25) gives completely inaccurate predic- 
tions of C p . However, the simple wall function with 
n -1 " < 50 gives similar predictions compared with any 
of the other models operated at much smaller n+ val- 
ues (<2). Also, the simple wall function with a larger 
value of n + accelerated the convergence rate of the 
k-e model as shown in figures 17 and 18. These fig- 
ures show the L-2 norm of the residual RMS and of q 
as a function of the number of iterations, respectively. 
A solution with the wall function and n + = 50 was 
established in less than 2000 iterations, which was 
similar to the number of iterations required with the 
Baldwin-Lomax turbulence model. The following ta- 
ble summarizes the c^ and q predictions based on 
the different models and modifications. 


Procedure 

ci 

Cd 

Experiment 

0.743 

0.0242 

k-e 1, < 2 

0.720 

0.0257 

k-e 2, n + < 2 

0.764 

0.0269 

fc-£3, n + <2 

0.772 

0.0257 

Baldwin-Lomax, n + < 2 . . . 

0.756 

0.0290 

WF, n + < 25 

0.730 

0.0261 

WF, n+ < 50 

0.736 

0.0231 

k-el, n + < 25 

0.364 

0.0222 


Supersonic Flow Through Square Duct 

Numerical calculations that used linear and non- 
linear k-e turbulence models were carried out for 
supersonic flow through a square duct. Figure 19 
shows a schematic of the known secondary flow pat- 
tern in square duct flows. The flow is symmetric 
about the y- and 2 -axes so only one quadrant of the 
duct flow was computed. A Mach number of 3.9 
and a unit Reynolds number of 0.012 x 10 6 cm -1 
(0.035 x 10 6 in’” 1 ) were used with a 41 x 41 grid in the 
cross-flow plane and 251 grid points in the streamwise 
direction. Because the flow is complex, appropriate 
grid spacing near solid boundaries was maintained 
to ensure appropriate near-wall effects of k and e. 
The first point located off the wall was < 1 and 
the grid was stretched in the normal direction by an 
exponential grid-stretching formula. Approximately 
16 points were placed in each direction normal to the 
walls to resolve the boundary layer. The rest of the 
points in the normal direction were distributed uni- 
formly between the edge of the boundary layer and 
the symmetry boundary. 

Figures 20-24 show comparisons of the results ob- 
tained for the linear and nonlinear turbulence models 
and the experimental data. (See ref. 14.) Figure 20 
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presents the effects of linear and nonlinear turbu- 
lence models on the computed skin-friction distri- 
bution at x/D = 50. The skin- friction distribution 
obtained with the nonlinear model is in better agree- 
ment with the experimental data, whereas the linear 
model shows a monotonic increase in value for Cj 
with the spanwise coordinate. The nonlinear model 
captures the undulations observed in the experimen- 
tal data with remarkable precision. These undula- 
tions represent the convecting effect of the secondary 
flow, corrections which were not predicted by the lin- 
ear model. 

Figure 21 shows the effect of two wall-damping 
functions (eqs. (18) and (20)) on skin- friction pre- 
dictions when the nonlinear model is used. This re- 
sult demonstrates that to obtain improved predic- 
tions, the nonlinear model has to be coupled with 
the appropriate damping function such as that devel- 
oped by Speziale, Abid, and Anderson. (See ref. 17.) 
Previous studies with this near-wall model gave im- 
proved predictions in simple wall-bounded separated 
flows such as the backward- facing step. (See ref. 22.) 

Figure 22 shows the effect of including compress- 
ibility corrections (eqs. (16) and (17)) on skin- friction 
predictions with the nonlinear model. Compressibil- 
ity correction is clearly not needed for this case; the 
results obtained without a compressibility correction 
and with the Wilcox model (ref. 7) are both in close 
agreement with the experimental data. However, the 
Sarkar et ah model (ref. 6) gave lower values of Cf. 
These lower results are not surprising because the 
compressibility correction in reference 6 is applied 
to all regions of the flow and does not have a built- 
in mechanism to switch off near the wall boundaries 
where the compressibility effects are minimal. In 
contrast, the Wilcox compressibility correction has 
a switching function based on the local turbulence 
Mach number and turns off automatically in the re- 
gions with little or no compressibility. 

Figure 23 shows the effect of grid resolution on 
the computed skin-friction distribution. These com- 
putations were performed for 21, 31, and 41 points 
in the normal direction but with the same stream- 
wise grid spacing. When the mesh in the normal 
direction was refined, the stretching coefficient was 
progressively increased to obtain finer mesh spacing 
near the wall, but an identical number of grid points 
in the boundary layer was maintained. The points 
in the outer region increased by a factor of 3 in the 
transition from a coarse to a finer mesh. 

Figure 24 shows the cross-flow velocity patterns 
computed with the linear and nonlinear models and 
the experimentally measured cross-flow pattern at 


x/D = 50. Dramatically improved results are 
obtained with the nonlinear model shown in fig- 
ure 24(b). The results clearly show that the sec- 
ondary flows (vortices) are symmetrical about the 
diagonal and rotate in opposite directions. These 
vortices are essentially driven by the gradients of 
the Reynolds stresses, which cannot be simulated 
with the linear models and which transport net mo- 
mentum toward the corner of the duct. The com- 
puted cross-flow velocity vectors that are based on 
the nonlinear turbulence model agree well with the 
experimentally observed patterns. (See fig. 24(c).) 
In contrast, the linear model (fig. 24(a)), predicts 
a unidirectional flow because the turbulence model 
cannot adequately represent the flow physics. 

Concluding Remarks 

A systematic investigation was conducted to as- 
sess the effect of grid resolution and various near-wall 
damping functions in the turbulence model of kinetic 
energy (fc) and dissipation rate (s) on the computed 
flow field of several aerodynamic configurations. The 
computed results were compared with the available 
experimental databases. The geometries considered 
in the present study were a flat plate, 16° and 24° 
compression corners, a two-dimensional airfoil sec- 
tion, and supersonic flow through a square duct. In 
addition, a nonlinear k-e turbulence model was used 
to predict the aerodynamic characteristics of super- 
sonic flow through a square duct, and the effect of 
compressibility corrections was investigated. Skin- 
friction, pressure, and velocity distributions were rea- 
sonably predicted with the two-equation turbulence 
model in its different forms. 

The flat-plate test case was selected because it is 
the simplest of all the geometries for which the effects 
of various near- wall grid spacing, turbulence models, 
and damping functions can be easily tested. The 
results obtained with the Jones and Launder turbu- 
lence model {k-el) were more accurate and consistent 
than those for the Baldwin-Lomax turbulence model 
for large n + (height of first grid line normal to solid 
surface). The number of grid points normal to the 
solid surface is not as important as a correct value 
of n + . Use of the simple wall function with the k-e 
turbulence model allowed the use of coarser grids in 
which n + < 50 but which maintained reasonably ac- 
curate predictions of skin friction. 

The compression corner and airfoil geometries 
represent the next level of flow complexity because 
these cases contain separated flow regions that in- 
teract with a shock. A flow of Mach 2.85 over 16° 
and 24° compression ramps was computed and com- 
pared with the experimental data. For the case of 
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16° ramp, good agreement was obtained in predicting 
the surface pressure and the skin friction upstream 
of the ramp. The 24° ramp had a massive separated 
flow region. The k-e 1 model accurately predicted the 
reattachment point but predicted the location of the 
separation point slightly downstream of the experi- 
mental data. 

For the airfoil geometry, the k-e model with 
Jones and Launder damping function (k-e :1) and 
the Baldwin-Lomax turbulence model gave similar 
pressure coefficient distributions, although the k-e 1 
model predicted the shock slightly farther, in better 
agreement with the experiment. Both models also 
yielded a similar prediction for skin friction upstream 
of the shock. The simple wall function with the k-e 1 
turbulence model at large values of n + brought the 
local skin friction within the range of the predictions 
obtained for other models (predicted at < 2). 
Also, a larger n + accelerated the convergence rate 
of the k-e model. Comparisons were also made be- 
tween the models of Jones and Launder (fc-el), Van 
Driest (k-e2), and Speziale et al. (fc-^3). All models 
predicted similar pressure distributions, but use of 
different forms of the damping function yielded sig- 
nificantly different skin-friction predictions upstream 
of the shock. In general, the k-e 1 model most accu- 
rately predicted the overall features of the flow field. 

For the test case that featured square duct ge- 
ometry, a secondary flow structure developed in the 


direction perpendicular to the main flow. The skin- 
friction distribution with the nonlinear k-e 3 turbu- 
lence model was in better agreement with the ex- 
perimental data, whereas the linear k-e 3 turbulence 
model produced a monotonic increase in the local 
skin friction with the spanwise coordinate. The non- 
linear model clearly captured the major trends ob- 
served in the experimental data and flow-field fea- 
tures. The undulations that were observed represent 
the convecting effect of the secondary flow, undula- 
tions which were not predicted by the linear turbu- 
lence model. The near- wall damping function devel- 
oped by Speziale et al. (k-e 3) helped to yield better 
prediction of the skin-friction distribution than that 
of Jones and Launder (fc-el) for this case. The com- 
pressibility correction of Wilcox performed better 
than that of Sarkar et al. for this case because the 
former has a built-in mechanism to switch off near 
the wall where compressibility effects are small. 

This investigation provided significant insight into 
the applications of turbulence models in the predic- 
tion of attached and separated flows. Comparisons of 
grid effects and the use of different turbulence models 
indicate that the k-e turbulence model can be used 
successfully to predict these flows. 

NASA Langley Research Center 
Hampton, VA 23681-0001 
October 26, 1994 
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Appendix 


Special Cases 


Universal High-Order Symmetry 
Boundary Condition for Navier-Stokes 
Codes 


A universal symmetry boundary condition was 
derived by using a general three-dimensional curvi- 
linear coordinate system. Here, x\, x<i , and x 3 are 
the computational domain axes and U 2 , and u 3 
are the corresponding velocities. At each cell face, a 
local Cartesian coordinate system (Xu X 2 , and X 3 ) 
is established by the inward normal of that cell face 
and two other linearly independent directions in the 
plane of the cell face. A velocity vector in this neigh- 
borhood is decomposed into vector components in 
the local coordinate system. Across this boundary, a 
ghost point b is assumed to be located at the mirror 
image of the real point p within the computational 
domain. (See fig. A(l).) The velocity vector asso 
ciated with the ghost point is the mirror image of 
the original velocity vector and satisfies the follow- 
ing condition: 



(Al) 


Halfplane Symmetry 

In the case of halfplane symmetry, the symmetry 
plane is aligned with the xi~x 3 plane and the N array 
is defined as 

N = {0,l,0} r 

Substituting the values for N in equations (A2) 
and (A4) yields 


pi] 

l & 1 

f 1 

u 2 , 

H 

1 -U2 \ 

l u 3 J 

1 1 

l u 3 j 


Quarterplane Symmetry 


In the case of quarterplane symmetry, two sym- 
metry planes exist. One is the same as the halfplane 
symmetry case; the other symmetry plane is aligned 
with the X 2 , 2:3 plane. The N array is defined as 

N = {l,0,0} r 


where Ui is related to the local velocities u t for points 
p and b through the A (transformation) matrix, 


U 6 = A u b 1 
U p = Au p J 


(A2) 


where 



' n\ 

ri2 

^3 


A = 

m\ 

m 2 

7713 


A -1 = . 

_ *1 

A t 

h 

h m 



(A3) 


and n, m, and l are the directional cosines between 
the X and x coordinates. Equation (A2) can be 
rewritten from equations (Al) and (A3) as follows: 


u b = A~ l V b 


= A- 1 (U P -2U 1 {1,0,0} T ) 

= A _1 Au p - 2UiA r {l,0,0} r 


Furthermore, the above equation can be written in 
the simpler form, 

u b — u p — 2Ui{ni,n2,n 3 } r (A4) 


Substituting the values for N in equations (A2) 
and (A4) yields 


{ m ( u\ y 

srU/ 

This mathematical formulation is implemented in 
the PAB3D Navier-Stokes code and replaces all the 
special cases of symmetry boundary conditions. It 
is also used for specifying a slip boundary condition 
for Euler flow calculations. The halfplane symmetry 
for a polar grid is an exception to this rule. The 
velocity vector image is reflected across the halfplane 
boundary, whereas the ghost points at the pole are 
reflected across a moving mirror. Nevertheless, the 
generalized symmetry boundary condition can be 
easily modified to work in this case, and it has 
been included in the PAB3D code. This generalized 
symmetry condition has tremendously simplified the 
boundary condition procedure in advanced Navier- 
Stokes and Euler CFD methods. It also provides the 
ability to include symmetry boundary conditions in 
a structured CFD grid without the restrictions for 
alignment with the global coordinate system. 
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Figure Al. Symmetry boundary condition formulated in local Cartesian coordinate system. 
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4 e = P k = f(e,k) 

Figure 1. Inflow boundary condition for k-e turbulence model. 
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Figure 2. Solution for subsonic and supersonic flows with k-e 1 over flat-plate boundary layer. 



Figure 3. Solution for supersonic flow over flat-plate boundary layer based on k-e 1 and Baldwin-Lomax models. 
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(b) Baldwin-Lomax turbulence model. 


Figure 5. Effect of spacing first grid point off wall on velocity profile. 




(b) Baldwin-Lomax turbulence model. 

Figure 6. Effect of spacing first grid point off wall on average skin-friction distribution 
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(b) Average skin- friction distribution. 

Figure 9. Solution for k-s model with and without wall functions for n 
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(b) Skin- friction distribution. 

Figure 10. Computations for 16° compression ramp with Baldwin-Lomax and k-e 1 models. 









(c) Velocity profile. 
Figure 12. Concluded. 



(a) Pressure distribution. 



(b) Local skin-friction distribution. 

Figure 13. Computations for two-dimensional airfoil section with Baldwin-Lomax and k-e 1 turbulence models. 
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Figure 15. Local skin- friction distribution for two-dimensional airfoil section for different forms of k-£ turbulence 
models with and without wall function. 



Figure 16. Wall pressure distribution for two-dimensional airfoil section for different forms of k-e turbulence 
models with and without wall f un ction. 
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Figure 17. Convergence history for k-el with and without wall function and Baldwin-Lomax turbulence models. 



Figure 18. Lift coefficient for fc-el with and without wall function and Baldwin-Lomax turbulence models. 
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Figure 20. Local skin-friction distribution with linear and nonlinear models for supersonic flow through square 
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Figure 21. Effect of near- wall terms on nonlinear model for supersonic flow through square duct. 





Figure 22. Effect of compressibility correction on nonlinear model for supersonic flow through square duct. 



Figure 23. Effect of grid resolution and k-e 3 model on nonlinear model for supersonic flow through square duct. 
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